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ABSTRACT 



A robust numerical method called the Preconditioned Bi-Conjugate Gradient 
(Pre-BiCG) method is proposed for the solution of radiative transfer equation 
in spherical geometry. A variant of this method called Stabihzed Preconditioned 
Bi-Conjugate Gradient (Pre-BiCG-STAB) is also presented. These are iterative 
methods based on the construction of a set of bi-orthogonal vectors. The ap- 
phcation of Pre-BiCG method in some benchmark tests show that the method 
is quite versatile, and can handle hard problems that may arise in astrophysical 
radiative transfer theory. 

Subject headings: Line: formation - radiative transfer - scattering - methods: 
numerical 
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Introduction 



The solution of transfer equation in 
even 75 years after the first attempts by 



spherical geom e try re n aains a c 



Chandrasekharl (jl934j ): 



assic problem 



KosirevI (119341 ). w ho used 



Eddington approxim ation. In later decades more accurate methods we r e given (see 



978 



1973 



Peraiah 



2002 



for historical reviews). 



Hummer fc Rybickil (jl97ll ) 



Mihalas 



Kunasz fc Hummer 



1974 ) developed a variable Eddington factor method, and computed the solution 
on rays of constant impact parameter (tangents to the discrete shells and parallel to line 
of sight) in ID spherical geometry. This is a very efficient differential equation based 
technique which uses Feautrier solution alor ig rays of cons t ant im pact parameter. An 



integral equation method was develop ed by 



again based on tangent rays approach. 



Peraiah fc Grant 



Schmid-Burgkl (119741) to solve the problem, 



(119731 ) presented a highly accurate 



finite difference method based on the first order form of the transfer equation. All these 
methods were later extended to expanding, and highly extended atmospheres. However, in 
this paper we confine our attention to static, ID spherical atmospheres. 

In a next epoch in the development of spherical radiative transfer, the integral operator 
techniques were proposed. The idea of operator splitting and the use of approximate 



operators in iterat i ve me th ods was b r ought to the astrophysical radiative transfer in planar 



media by 



Cannon! (119731 ). IScharmerl (jl98ll ) extended his work with a new definition of 



the approximate operator. The applicatio n of In tegra. 



transfer started with the work of 



Hamann 



(11985h and 



operator technique to t 



Werner fc Husfeldl (119851 ). They used 



le spherical 



approximate operators that are diagonal, constructed from core saturati on approach. The A 



operator contains the non-local coupling between all the spatial points. lOlson et al.l (119861 ) 
showed that the diagonal part (local coupling) of the actual A operator itself is an optimum 
choice for the 'approximate operator'. These methods are known as approximate Lambda 



Iteration (ALI) methods. The ALI methods which are based on the concept of operator 
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splitting and the use of Jaco bi iterative technique, were w idely used in the later decades in 



radiative transfer theory (see 



Hubeny 


2003; 


Hamann 


2003 



for historical reviews). 



Gros et al. 



(119971 ) used an implicit integral method to solve static spherical line transfer 



problems. The most recent and interesting work on sp herical radiative tran s 



papers by 



Asensio Ramos fc Trujillo Buend (120061 ) and 



Daniel fc Cernicharol t00& \ both 



'er are the 



of which are based on Gauss-Seidel (GS) and Successive Over Relaxation (SOR) iterative 
techniques. 



Klein et al.l (119891 ) were the first to use BiCG technique in astrophysics. They use BiCG 
with incomplete LU decomposition technique in their double splitting iterative scheme along 
with Orthom in acceleration. They applied it to mult i- dimensional line transfer problem. 



Auerl (119841 ) describes a variant of Orthomin acceleration which uses 'Minimization with 
respect to a set of Conjugate vectors'. He uses a set of n (usually n=2 or n=4) conjugate 
direction vectors which are orthogonal to each other, constructed using the residual vectors 
with a purpose to accelerate the convergence sequence. 



Hubeny fc Burrows 



(120071 ) developed GMRES (actually its variant called Generalized 
Conjugate Residuals GCR) method to solve the spherical transfer problem. It is based on 
an application of the idea of Krylov subspace techniques. They applied it to a more general 
time-dependent transport with velocity fields in a medium which scatters anisotropically. 
They apply GMRES method to the neurino transfer. It can also be used for radiation 
transfer problem, including the simple problem of 2-level atom line transfer discussed in 
this paper. 



The Preconditioned Bi-Conjugate Gradient method (hereafter 



2000r) was first introduced to the line transfer in planar media, by 



're-BiCG, see eg.. 



Saad 



Paletou fc Anterrieu 



(12009! ) who describe the method and compare it with other prevalent iterative methods. 



namely GS/SOR. In this paper we adopt the Pre-BiCG method to the case of spherical 
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media. We also show that the 'Stabihzed Preconditioned Bi-Conjugate Gradient (Pre- 
BiCG-STAB)' is even more advantageous in terms of memory requirements but with similar 
convergence rate as Pre-BiCG method. 

It is well known that the spherical radiative transfer in highly extended systems, 
despite being a straight forward problem, has two inherent numerical difficulties namely (i) 
peaking of the radiation field towards the radial direction, and (ii) the (1/t^) dilution of 
radiation in spherical geometry. To handle these, it becomes essential to take a very large 
number of angle {^) points and spatial (r) points respectively. The existing ALI methods 
clearly slow down when extreme angular and spatial resolutions are demanded (for example 
see table [1]). Therefore there is a need to look for a method that is as efficient as ALI 
methods, but faster, and is relatively less sensitive to the grid resolution. The Pre-BiCG 
method and a variant of it provide such an alternative as we show in this paper. 

Governing equations are presented in Sect. 12. 1[ In Sect. \2.2\ we define the geometry 
of the problem and the specific details of griding. In Sect. 12.31 the benchmark models 
are defined. We briefiy recall the Jacobi, and GS/SOR methods in Sect. 12.41 In Sect. [3] 
we describe the Pre-BiCG method. The computing algorithm is presented in Sect. 13.11 
In Sect. IHwe describe the Pre-BiCG-STAB method briefiy, and we give the computing 
algorithm in Sect. 14.11 In Sect. [5] we compare the performance of Pre-BiCG with the 
Jacobi, and GS/SOR methods. In Sect. [6] we validate this new method, by comparing with 
the existing well known benchmark solutions in spherical line radiative transfer theory. 
Conclusions are presented in Sect. [71 
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2. Radiative transfer in a spherical medium 

2.1. The transfer equation 

In this paper we restrict ourselves to the case of a 2-level atom model. Further, we 
assume complete frequency redistribution (CRD). The transfer equation in divergence form 
is written as 

dlir, fi,x) 1 — /i^ dl(r, /i, x) 
II- 



dr r djj, 

= [XLir)(t){x) + Xcir)] 

x[S{x,r) - I{r,n,x)]. (1) 

Here, / is the specific intensity of radiation, S - the source function, r - the radial distance, 
fi - the direction cosine, x - the frequency measured in Doppler width units from line center, 
0(x) - the line profile function, and Xl(^), Xc{r) - line center and continuum opacities 
respectively. The differential optical depth element is given by 

dr(r) = -XLir)dr. (2) 



Ther e are several methods which use the above form of the transfer equation (see iPeraiah 



20021 ). In our paper we solve the transfer equation on a set of rays tangent to the spherical 
shells. It is written as 

±^^^|^ = [XL(r)0(x) + Xc(r)] 

x[S{x,r)-I^{z,p,x)], (3) 

for the outgoing (+) and incoming (-) rays respectively. Here z is the distance along 
the tangent rays and p is the distance from the center to the points on the vertical axis 
(the mid-line), where the tangent rays intersect it (see Fig. [T]). The direction cosines 



(0 < /i < 1) are related to p by /i = \/l — (p^/r^) for a shell of radius r. The optical 
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depth scale along the tangent rays are now computed using dr(z) = dr(r) / yU. In practical 
work, due to the symmetry of the problem, it is sufficient to perform the computations on 
a quadrant only. The source function is defined as 

cf X XL(r)0(x)5L(r) + Xc(r)5c(r) 

Sc{r) is the continuum source function taken as the Planck function By{r) throughout this 
paper. The monochromatic optical depth scale Ar^, = Arz[0(x) +/?c], with jSc = Xc{^) / Xhi^) 
along the tangent rays. For simplicity, hereafter we omit the subscript z from and write 
r to denote t^. The line source function is given by 

S,(r) = (l-e) / '^ r dx' 



-1 2 

0(x')/(r,/i',x') + eB,(r), (5) 

with the thermalization parameter defined in the conventional manner as e = Cui/ (Aui + Cui), 
where Cui and are collisional and radiative de-excitation rates. The intensity along the 
rays is computed using the formal solution integral 

I^{t,p,x) = /o^(r,p,x)exp[-Ar^] + 

^ exp[-Ar^]5(r')[0(a;)+/3c]dr'. (6) 
The corresponding integral for the incoming rays is 

I'{t,p,x) = /o"(r,p,a;)exp[-Ar^] + 

r exp[-Ar'JS{T'Mx)+[3,]dT'. (7) 
Jo 

Here, Iq{t,p,x) represents the inner boundary condition imposed at the core and along the 
mid- vertical line (see Fig. [1]). Iq{t,p,x) is the outer boundary condition specified at the 
surface of the spherical atmosphere. When the above formal integral is applied to a stencil 
of short characteristic (MOP) along a tangent ray, it takes a simple algebraic form 

I^{t,p,x) = /^(r,p,a;)exp[-ArM] + 
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*±(r,p,x)5p, (8) 

where S'm.o.p are the source function values at M, O and P points o n a short charac t eristic . 
The coefficients ^ are calculated following the method described in lKunasz fc Auerl ( 119881 ) . 



2.2. The constant impact parameter approach 

In Fig. [H we show the geometry used for computing the specific intensity I{t,p,x) 
along rays of constant impact parameter. 

In a spherically symmetric medium, we ffist discretise the radial co-ordinate r 
(-Rcore r < R), whcrc -Rcore IS the core radius, and R is the outer radius of the atmosphere. 
The radial grid is given by rk, = 1, 2, . . . , A^d, where ri is the radius of the outer most 
shell, and is that of the inner most shell. dQ/An is the probability that the direction 
of propagation of an emitted photon lies within an element of solid angle dQ. In the 
azimuthally symmetric case, it is dfi/2. To calculate the mean intensity J in plane parallel 
geometry, we integrate the intensity over the angular variable fi itself. In spherical medium, 
we have one to one correspondence between the (/i, r) and the {p, r) system. In {p, r) 
system, the probability that a photon is emitted with its impact parameter between p and 
p+dp, propagating in either positive /i or negative /i direction is pdp/2r\Jr'^ — p"^. The 
direction cosines made by the rays in the (yU, r) space, with a tangent ray of constant p 



value, are given by /i; = a/1 — (p'^/rf) at different ra dii r-,. Therefore the 



factor dn/2 can be changed to pdp/2r^^/r'^ — p"^ (see lKunasz fc Hummeiill973l ). 



angu 



ar integration 



The p - grid construction: If A'^c is the number of core rays, then the p-grid for the core 
rays is computed using: 



do i = 1, A'c 
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< P < Rcove 

end do 

The number of lobe rays equals the number of radial points. For lobe rays, the p-grid is 
same as radial r-grid. It is constructed using: 

do z = 1, A^d 
p{Nc + i) = r{i) 
end do 

where A'^ti=the number of radial points . Thus , the total number of impact parameters is 



Np = Nc + Na- We have followed 



Aued (119841 ) in defining the p-grid in this manner. 



2.3. Benchmark models 

Geometrical distances along the rays of constant impact parameter are constructed as 
follows: 

z{p,r) = A/r^ — j9^. (9) 

For spherical shells we perform several tests using power-law type variation of density. For 
such atmospheres, the line and continuum opacities also vary as a power law given by 

XL,c(r) (xr-". (10) 

Let C and C denote the proportionality constants for Xhi^) and Xc{^) respectively. The 
constant C can be determined using the optical depth at line center T. For a power law 
with index h, 



C= ^<'-»),, (11) 
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using the given input value of /3c = C/C we can compute the constant C. 

We use Voigt profile with damping parameter a or the Doppler profile for the results 
presented in this paper. The spherical shell atmosphere is characterized by the following 
parameters: (i?, n, T, a, e, (3c-, B^). We recall that R is the outer radius of the spherical 
atmosphere surrounding a hollow central cavity of radius -Rcore- When R = -Rcore we recover 
the plane parallel limit. For the spherical shell atmospheres, we take -Rcore = 1 as the unit 
of length to express the radial co-ordinate. The boundary conditions are specified at the 
outer boundary (/^(t = 0,p,x) = 0) and the inner boundary. There are two types of inner 
boundary conditions: 

(a) Emitting Core: 
Core rays: 

I^iT = T,p< Rcore, X)=B,. (12) 

Lobe rays: 

/+(r = T,p = ri,x) = /"(r = T,p = ri,x), (13) 
i = 1, 2, . . . , iVd along the mid- vertical. 



(b) Hollow Core: 

For both the core and the lobe rays: 



I^{t = T,p,x) = I {t = T,p,x). 



(14) 



The hollow core b oundary condition is also called 'planetary nebula boundary condition' 



see 



Mihalas 



19781 ). It is clear that a spherical shell with a hollow core is equivalent to a 
plane parallel slab of optical thickness 2T with symmetry about the mid-plane at r = T. 
We use spherical shell atmospheres for most of our studies. 
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2.4. Iterative methods of ALI type for a spherical medium 



The ALI methods have been su ccessfully used for the solution of transfer equation in 



spherical shell atmospheres (see eg., 



Hamann 



20031 . a nd references therein). These authors 



use the Jacobi iterative methods (first introduced by 



Olson et al 



(Il986l )) for computing 



the source function cor rections. Recently the GS method has 



same problem (see e g., 



Hubeny fc Burrows 



Asensio Ramos fc Trujillo Bueno 



2006 



jeen proposed to solye the 



Daniel &: Cernicharo 



20081 ). 



(120071 ) proposed the GMRES method for solving spherical radiative 
transfer problem. GMRES and Pre-BiCG both belong to Krylov subspace technique. In 
this paper we compute the spherical transfer solutions by Jacobi and GS/SOR methods, 
and compare with the solutions computed using the Pre-BiCG method. For the sake of 
clarity, we recall briefiy the steps of Jacobi and GS/SOR methods. 



Jacobi Iteration Cycle: The source function corrections are given by 

Xan _ an+l an _ ~ ^)'^k + ~ '^k /i p-N 

[1 - (1 - ^)K,k] 

for the nth iterate. Here k is the depth index. The A* is the approximate operator which is 
simply taken as the diagonal of the actual A operator defined through 



A[S] = J; (16) 

^ d^' 
1 ^ 



J{t)= I ^/ dx'(f){x')I{T,i2\x'). (17) 



oo 



GS/SOR Iteration Cycle: The essential difference between the Jacobi and GS/SOR methods 
is the following: 

5^+1 = + (18) 
Here the parameter u is called the relaxation parameter which is unity for the GS technique. 
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The SOR method is derived from th e GS method by simply taking 1 < < 2 (see 



Trujillo Bueno fc Fabiani Bendicho 



the GS method is given by 



SSI 



19951 . for details). The source function correction for 

~n(old+new) j-, ryn 



1 - e)J. 



(19) 



[1 - (1 - e)As:,k] 

where the quantity denotes the mean intensity computed using new values of 

the source function as soon as they become available. For those depth points for which the 
source function correction is not yet complete, GS m ethod uses the values of the source 



function corresponding to the previous iteration (see 



Trujillo Bueno &: Fabiani Bendicho 



19951 ). For clarity we explain how the GS algorithm works in spherical geometry, on rays of 



constant impact parameter. 



Begin loop over iterations 

Begin loop over radial shells with index k 

Begin loop over impact parameters (or directions) with increasing p 



For the nth iteration: 

For the incoming rays (yU < 0): 
(Reverse sweep along radial shells) 

(a) This part of the calculations start at the outer boundary for all impact parameter rays. 



(b) Jk are first calculated for a given radial shell k using S^, S^_^ and S^^ 



(c) The partial integral Jk(/i < 0) are calculated before proceeding to the next shell. This 
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part of the calculations is stopped when the core (for the core rays) and the mid- vertical 
line (for the lobe rays) are reached. 

For outgoing rays (/x > 0): 
(Forweird sweep along radial shells) 

(d) This part of the calculations start at the inner boundary. First, for the radial shell with 

k = Nd 

Jnjj is calculated, using boundary conditions /Na- 

(e) SS^^ is computed and the source function is updated using S^'^^ = S^^ + SS^^. 

(f) For the next radial shell k = — 1: 

to calculate I-^^-i by applying the short characteristic formula, Snj, 'S'Nd-i ^-nd <S^Nd-2 
are needed. Already S^'^^, "^Na-i available. GS takes advantage of the 

available new source function ai k = N^^. I^^-i is calculated with this set of source functions. 

(g) Then JNj_i(/i > 0) are calculated using -fNd-i- 

(h) Note that, JNd-i(Ai < 0) was calculated using S^^, S^^_i, and S^^_2 whereas 
JNd-i(A* > 0) used the "updated" source function S^~^^. Therefore JNd-i is corrected by 
adding the following correction: 

AJNd-i = SS^^ J'^ *^,(// < 0) d//. 

(i) (^-S^Nd-i "^n"*^-! — "^Nd-i + <^'^Nd-i calculated. 
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(j) Since "updated" ai k = — 1 is also available now, before going to the next 

radial shell it is appropriate to correct the intensity at the present radial shell by adding to 
it, the following correction term 

End loop over impact parameters (or directions) 
End loop over radial shells 
End loop over iterations. 



3. Preconditioned BiCG method for a spherical medium 



f the Pre-BiCG method. The 



In this section we first describe the essential ideas o: 
complete theory of the method is described in 
source function with a background continuum is given by 

XL(r)0(x)5L(r) +Xcir)Scir) 



SaadI fl2000h . We recall that the 2- level atom 



S(x, r) 



XL(r)0(x) + xc(^) 



It can be re-written as 



(20) 



S{x, r) = p{x, r)Si,{r) + (1 — p{x, r))Sc{r) 



where 



p{x, r) 



Xh{r)(t){x) 



X-Lir)(t){x) + xcir) 



From equations ([5]), (jT6]l and (jlTll . we get 



5(x, r) = p(x, r){(l - e)A[5(x, r)] + eB^(r)} + 
(1 -p(x,r))5c(r). 



(21) 



(22) 



(23) 
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Therefore the system of equations to be solved becomes 

[/ — (1 — e)p{x, r)A]S{x, r) = p{x, r)eBy{r) + 

(l-p(x,r))5e(r), (24) 

which can be expressed in a symbohc form as 

Ay = b; with 

i = [/ - (1 - e)p{x, r)A]; y = S{x, r). (25) 

The vector b represents quantities on the RHS of equation fl2^ . Now we describe briefly, 
how the Pre-BiCG method differs from ALI based methods. 

Let denote the n-dimensional Euclidean space of real numbers. 

Definition: The Pre-BiCG algorithm is a process involving projections onto the 
m-dimensional subspace (m < n) of M" 

Kra = sp8in{vi,Avi,...,A'^-^Vi}, (26) 

and also being orthogonal to another m-dimensional subspace of M" 

£m = span{iOi, A^wi, A^^'^-^^Wi}. (27) 

Here Vi is taken as the initial residual vector Vq = b — Ayo with yo the initial guess for 
the solution of equation fl2^ . The vector Wi is taken as arbitrary such that the inner 
product {vi,Wi) 7^ 0. The method recursively constructs a pair of bi-orthogonal bases 
{vi] i = 1,2 . . . , m} and {wi, i = 1,2..., m} for /Cm and C^n respectively, such that they 
satisfy the bi-orthogonality condition {vi,Wj) = 6iy For the purpose of application to the 
radiative transfer theory it is convenient to write the Pre-BiCG steps in the form of an 
algorithm For simplicity we drop the explicit dependence on variables. 
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3.1. The Preconditioned BiCG Algorithm 



Our goal is to solve equation (125|) . In this section the symbols Vi and Pi are used to be 
in conformity with the standard notation of residual and conjugate direction vectors. They 
should not be confused with the radius vector ri and impact parameter pi which appear in 
spherical radiative transfer theory. 

(a) The very first step is to construct and store the matrix (which does not change with 
iterations, for the cases considered here, namely 2- level atom model). Details of computing 
efficiently is described in appendix Rl 

We follow the preconditioned version of the BiCG method. Preconditioning is a 
process in which the original system of equations is transformed into a new system, which 
has faster rate of convergence. For example, this can be done by solving the new system 
M-^Ay = M~^b where M i s an a ppropriately chosen matrix, called the "preconditioner" 



(See also equation 2 of 



Auer 



19841 ). This preconditioner is chosen in such a way that, 

(i) the new system should be easier to solve, 

(ii) itself should be inexpensive to operate on an arbitrary vector, 

(iii) the preconditioning is expected to increase the convergence rate. 

The choice of the preconditioner depends on the problem at hand. When an appropriate A* 
is chosen such that t he amplification rn atrix [/ — (1 — e)A*] has as small a maximum eigen 



value as possible (see 



Olson et al 



19861 ). the convergence rate is enhanced. What enables 
the convergence of ALI, that satisfies the above property, and simplest to manipulate, is 
the diagonal of the A itself. Therefore the amplification matrix [/ — (1 — e)A*] with a 
diagonal form for A* is a simple and natural choice as a 'preconditioner'. We construct the 
preconditioner matrix M by taking it as the diagonal of A. 



(b) An initial guess for the source function is 



(28) 



where the thermal part eB is taken as an initial guess for Sl- 



(c) The formal solver is used with Uq as input to calculate J{yo)- 



(d) The initial residual vector is computed using 



ro = 6 - Ayo. 



(29) 



(e) The initial bi-orthogonal counterpart for rg is chosen such that we have (ro, rg) 7^ 0. 
One can choose Tq = Vq itself. 

Such an initial choice of Vq vector is necessary, as the method is based on the 
construction of bi-orthogonal residual vectors ri and recursively, for i = 1,2, ... ,m, 
where m is the number of iterations required for convergence. The process of constructing 
the bi-orthogonal vectors gets completed, once we reach the convergence. In other words, 
the number of bi-orthogonal vectors necessary to guarantee a converged solution represents 
the actual number of iterations itself. It is useful to remember that when we refer to 
'bi-orthogonality' hereafter, say eg., of the residual vectors Vi, r* we simply mean that 
(ri, fj*) — ioT i j, but (ri, r*) need not be unity. 

(f) The bi-orthogonalization process makes use of conjugate direction vectors p 
and p* for each iteration. They can be constructed during the iterative process, again 
through recursive relations. An initial guess to these vectors is made as po — Vq and Pq — Vq. 
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(g) The preconditioned initial residual vectors Co computed using 

Co* = M-V*. (30) 

(h) For i — 1,2, ... , the following steps are carried out until convergence: 

(i) Using the formal solver with pi as input (instead of actual source vector y), J[pi\ is 
obtained. 

(j) A\pi] is computed using 

i[pi]=Pi-(l-e)pJ[pi]. (31) 

(k) The inner products 

{A[p,],p*) and {vuCn, (32) 
are computed and used to estimate the quantity 

(1) The new source function is obtained through 

= 2/i + oiiPi. (34) 

Test for Convergence: Let uj denote the convergence criteria. If 

max{6y/y} < uj, (35) 

T 

then iteration sequence is terminated. Otherwise it is continued from step (m) onwards. 
The convergence criteria lj is chosen depending on the problem. 
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(m) Following recursive relations are used to compute the new set of vectors to be used in 
the {i + l)th iteration: 

rj+i = r; - (36) 
r*^, = r*-a,A'[p*], (37) 
C+i = M-'rl,. (38) 

(n) The quantity /5i is computed using 

A ^ (39) 

(o) The conjugate direction vectors for the {i + l)th iteration are computed through 

Pi+i = "Ti+i + PiPi, 

K+i = ^iVi + Apr- (40) 

(p) The control is transferred to step (g). 

The converged source function y is finally used to compute the specific intensity 
everywhere within the spherical medium. 



4. Transpose free variant - Pre-BiCG-STAB 

In spite of higher convergence rate, computation and storage of the matrix is a 
main dis-advantage of the Pre-BiCG method. To avoid this, and to make use of only the 
'acti on' oi A ma trix on an arbitrary vector, a method called 'BiCG-squared' was developed 



(See ISaadI l200d . for references and details), which is based on squaring the residual 
polynomials. Later it was improved by re-defining the residual polynomial as a product of 
two polynomials and obtaining a recursive relation for the new residual polynomial. This 
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product involves residual polynomial of the Pre-BiCG method and a new polynomial which 
'smoothens' the iterative process. In this section we give the computing algorithm of the 
Pre-BiCG-STAB method as applied to a radiative transfer problem. As described below, 
we can avoid computing and storing of the matrix in the Pre-BiCG-STAB method. 
However we would now need to call the formal solver twice per iteration unlike in Pre-BiCG 
method, where it is called only once. This results in an increase in number of operations 
per iteration when compared to Pre-BiCG method, causing a slight increase in the CPU 
time per iteration. In spite of these the Pre-BiCG-STAB method turns out to be always 
faster than the regular Pre-BiCG method in terms of convergence rate (lesser number of 
iterations for convergence). 

4.1. Pre-BiCG-STAB algorithm 

Now we give the algorithm of Pre-BiCG-STAB method to solve the system 
M~^AS — M~^b. Here M is a suitably chosen preconditioner matrix. The computing 
algorithm is organized as follows: 

(a) First initial preconditioned residual vectors and conjugate direction vectors are defined 
through 

Zo = - M-^AS, (41) 

z* = Zo, Po = Zq. (42) 

(b) For j — 1,2,... the following steps are carried out until convergence. 

(c) Using Pj instead of the source function a call to the formal solver is made to compute AP^. 



(d) The coefficient ctj can be evaluated now as 

^ ^ ^^j'/o) (43) 

(e) Another vector qfj is calculated as 

= Zj - ajM-^Pj. (44) 

(f) Using Qj in place of the source function a call to the formal solver is made to obtain Aqy 

(g) The coefficient cuj is estimated as 

. .<^:'^^i-»i>, . (45) 

(h) The updated new source function is calculated as 

5j+i = 5j + ajPj+Wjqj. (46) 

(i) Test for convergence is made as in the Pre-BiCG algorithm. 

(j) Before going to the next iteration a set of recursive relations are used to compute 
residual vectors 

Zj+i = gj-WjM-^iqj, (47) 

and conjugate direction vectors 

Pj+i = Zj+i + /3j(Pj - ^jM-^iPj). (48) 
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for the next iteration, where the coefficient /?j is 



(k) The control is now transferred to the step (b) . 



5. Compcirison of ALI and Pre-BiCG methods 

There are two characteristic quantities that define iterative techniques. They are (a) 
convergence rate, which is nothing but the maximum relative change (MRC) defined as 

it:e = max{^}, (50) 

and (b) the total CPU time Ttotai required for convergence. Ttotai is the time taken to reach 
a given level of convergence, taking account only of the arithmetic manipulations within 
the iteration cycle. We also define a quantity called the true error and use it to evaluate 
these methods. 



5.1. The behaviour of the mciximum relative change (MRC) 

In this section we compare Rc and Ttotai for the Jacobi, GS, SOR, Pre-BiCG and the 
Pre-BiCG-STAB methods. The SOR parameter used is 1.5. It is worth noting that the 
overrates (the time taken to prepare the necessary set up, before initiating the iterative 
cycle) are expected to be different for different methods. For instance, in Jacobi and 
GS/SOR this is essentially the CPU time required to set up the A* matrix. In the Pre-BiCG 
method this involves the time taken to construct the matrix, which is a critical quantity 
of this method. The Pre-BiCG method is described in this paper in the context of a 2-level 
atom model, because of which, we do not need to update the matrix at each iteration. 



- 23 - 



For the Pre-BiCG-STAB method it is the time taken to construct the preconditioner matrix 
M. 

Fig. [2] shows a plot of Rr for differe nt methods. We can take -Rc as a measure of the 



convergence rate. 



Chevalher et aL 



(120031 ) show that it always becomes necessary to use high 
resolution grids, to achieve high accuracy of the solution (See also Sect. [T]of this paper). 
This is especially true in the case of spherical radiative transfer where a spatial grid with 
a large number of points per decade becomes necessary to achieve reasonable accuracy. In 
the following we discuss how different methods respond to the grid refinement. It is a well 
known fact with the ALI methods, that the convergence rate is small when the resolution of 
the depth grid is very high. In contrast they have a high convergence rate in low resolution 
grids. On the other hand the Rc of Pre-BiCG and Pre-BiCG-STAB methods have higher 
convergence rate even in a high resolution grid. Fig. [2](a) shows Rc for different methods 
when a low resolution spatial grid is used (5pts/D in the logarithmic scale for r grid). The 
Jacobi method has a low convergence rate. In comparison, GS has a convergence rate 
which is twice that of Jacobi. SOR has a rate that is even better than that of GS. However 
Pre-BiCG and the Pre-BiCG-STAB methods have the higher convergence rate. Fig. [2](b) 
and Ws) are shown for intermediate (8 pts/D) and high (30 pts/D) grid resolutions. 
The essential point to note is that, as the grid resolution increases, the convergence rate 
decreases drastically and monotonically for the Jacobi and the GS methods. It is not so 
drastic for the SOR method which shows non-monotonic dependence on grid resolution. 
The Pre-BiCG and Pre-BiCG-STAB methods exhibit again a monotonic behaviour apart 
from being relatively less sensitive to the grid resolution. 

In Table [1] we show what happens when we set convergence criteria to progressively 
smaller values {uj =10~®, 10"*^, and 10~^° for Tables 1(a), 1(b) and 1(c) respectively) for 



various grid resolutions. The model used to compute these results is (n, i?, T, a, e, (3c, By) 
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(0, 10, 10^, 10-^ 10-^ 0, 1). The idea is to demonstrate that for a given grid resolution 



(corresponding rows of the Tables 1(a), 1(b) and 1(c)), all the methods show a monotonic 



increase in the number of iterations for convergence, as we decrease the uj. On the other 
hand Pre-BiCG and Pre-BiCG-STAB require much less number of iterations to reach the 
same level of accuracy. 

CPU time considerations: Table [2] shows the CPU time requirements for the methods 
discussed in this paper. The model used to compute these test cases is (fi, R, T, e, Pc, 
B^, u)= (0, 300, 10^, 10-^ 0, 1, 10-^). The grid resolution considered is 30 pts/D. The 
CPU time for convergence can be defined as the computing time required to complete 
the convergence cycle and reach a fixed level of accuracy. We recall that the overrates 
in computing time is the time taken to prepare the necessary set up before initiating the 
iterative cycle. Total computing time is the sum of these two. In appendix |X] we discuss in 
detail how to construct A* matrix for Jacobi, GS/SOR methods and A and M matrices for 
Pre-BiCG and Pre-BiCG-STAB methods respectively with an optimum effort. Construction 
of these matrices constitutes the overrates in computing time of each method. The first 
row of Table [2] shows that Pre-BiCG is the fastest to complete the convergence cycle. The 
reason why Pre-BiCG-STAB takes slightly longer time than Pre-BiCG is explained at the 
end of Sect. HI 

The second row of Table [2] shows that all methods except Pre-BiCG take nearly 8 
seconds as overrates for the chosen model. Pre-BiCG takes additional 3-4 seconds as explicit 
integrals are performed for computing off-diagonal elements also (Unlike the other methods 
where such integrals are performed only for diagonal elements). 

The last row of Table [2] shows that in terms of total CPU time requirement, the other 
methods fall behind the Pre-BiCG and the Pre-BiCG-STAB. Pre-BiCG seems to be a bit 
faster compared to Pre-BiCG-STAB for the particular model chosen. However it is model 
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dependent. For instance, as the contribution towards overrates increases, Pre-BiCG-STAB 
clearly stands out as the fastest method of all, discussed in this paper. 



5.2. A study of the True Error 

We now study the true errors in these methods (see Figure [3]). The model parameters 
are (n, i?, T, e, Pc, By)= (0, 10, 10^, 10~^, 0, 1). A coherent scattering limit is used. To 
define a true error, we need a so called 'exact solution'. Except for highly idealized cases, 
exact solutions do not exist. For practical purposes, the exact solution can be defined as 
a solution obtained on a spatial grid of resolution that is three times larger than the grid 
resolution of the model that we are interested in. Also, we extend the iteration until Rc 
reaches an extremely small value of 10^^^. The source function compute d in this way can be 
called ^(oo, oo) (fully converged solution on an infinite resolution) (see 



Aueretal 



19941 ) 



The source function at the nth iterate is denoted by S". We define the true error as 



rr I ^cxact I / r- 1 \ 

ig = max I I, (51) 



following 



Trujillo Bueno fc Fabiani Bendichd (119951 ). In Fig. [3](a) we show Tg computed 



for the Pre-BiCG method using three grid resolutions, namely 10 pts/D, 14 pts/D and 20 
pts/D. The plateau of each curve represents the minimum value of the true error reached for 
a given grid resolution. We notice that as the resolution increases, Tg gradually decreases in 
magnitude as expected. In Fig. [31(b) we show Tg computed for the Pre-BiCG-STAB method. 
The model parameters are same as in Fig. [31(a). Clearly, Pre-BiCG-STAB shows a smooth 
decrement of true error compared to Pre-BiCG, because of the smoothing polynomial used 
to define the residual vectors. In Fig. [3t^c) we compare the decrement of true errors for 
different iterative methods. The grid resolution chosen is 14 pts/ D with other model 
parameters being same as in Figs. [3](a), (b). The decrease of the true errors follows the 
same pattern in all the iterative methods, although the number of iterations required for 
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Te to reach a constant value (plateau) depends on the method. To reach the same level of 
true error, the Pre-BiCG and Pre-BiCG-STAB methods require considerably less number 
of iterations, when compared to the other three. 



5.3. A theoretical upper bound on the number of iterations for convergence 

in the Pre-BiCG method 

Suppose that A is an x matrix. The solution to the problem Ay = 6 is a vector 
of length N^. In an A^^d-dimensional vector space Vn^ the maximum number of linearly 
independent vectors is N^. Hence, there can at the most be orthogonal vectors in Vn^. 
The Pre-BiCG method seeks a solution by constructing orthogonal vectors. We recall that 
the residual counterpart vectors {r]", rg, . . . , r|^} constructed during the iteration process 
are orthogonal to the initial residual vector vq. Thus, when we reach convergence after M 
iterations, we will have a set of M + 1 orthogonal vectors {vq, rjj', rg, . . . , tJ/}- From the 
arguments given above, it is clear that M + 1 < A^d, namely in the Pre-BiCG method, 'the 
convergence must be reached theoretically in at the most A"d steps (or ite rations)'. This sets 



Hestenes fc Stiefel 



an up per limit to the number of iterations to reach convergence (see also 
19521 ). For example when the dimensionality of a problem is high (very large value of A^d), 
the Pre-BiCG method ensures convergence in at the most Ad iterations. A theoretical upper 
bound on the number of iterations also exists for the Pre-BiCG-STAB method, whereas 
the other methods do not have such a theoretical upper bound. In practice we find that 
Pre-BiCG and Pre-BiCG-STAB methods actually require much less number of iterations 
than Ad, even when A'd is large. 
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Results and discussions 



The main purpose of this paper is to propose a new method to solve the line transfer 
problems in spherically symmetric media. In this section we show some illustrative examples 
in order to com pare with the famous bench marks for spherical transfer solutions presented 
in the papers by lKunasz fc Hummed (11974 ). In Fig. H] we show source functions for different 
test cases. Fig. |4](a) shows the source functions for e = 10~^ and e = 10~^ for R = 300 
and T = 10'^. Other model paramet ers are (n, 0r, B„) = ( 2, 0, 1). We use a Doppler 



profile to compare with the results of 



Kunasz &: Hummer! (jl974j ). Plane parallel result is 



also shown for comparison. When e = 10~^ we observe that the thermalization is reached 
at the thermalization length for the Doppler profile namely 1/e = 100. When e = 10~^ 
thermalization does not occur. Clearly the minimum value of the source function is eB,y. 
For the large values of i? = 300 and opacity index n = 2, as r increases the opacity decreases 
steadily and the source function indeed approaches this minimum value near the surface 



layers. For this case, the departure of the source function from planar 
the su rface. It can be shown (dashed line of Fig. Hl^b), see also Fig. 3 of 



imit is severe near 



Kunasz &: Hummer 



(jl974j )) that this departure is not so acute when n = 0, but is more acute when n = 3 
(dash triple-dotted line in Fig. IH^b)). In Fig. IH^b) we plot source function for the same 
model as Fig. Hl^a) but for various values of n. For negative h, the distinction between S{t) 
vs. T curves for different n is small. For positive n, the effects are relatively larger (see 
dot-dashed and long dashed curves in Fig. 111(b)). 

In Fig. 111(c), we show source function variation for a range of spherical extensions R. 
We have chosen an effectively optically thin model (T, e) = (10®, 10^^°) because in such a 
medium, thermalization effects do not completely dominate over the effects of sphericity. 
Other parameters are same as in Fig. Hl^b). Clearly, the decrease in the value of source 
function throughout the atmosphere is monotonic, with an increase in the value of R from 
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1 to 10^ 

In Fig. we show effects of limb darkening in splierical atmospfieres for R = 10^ and 
R = 10^. Tlie otlier model parameters for Fig. M^a.) are (n, T, e, /3c, B,^)= (2, 10^, 10"'^, 0, 
1). A Doppler line profile is used. From the Figure, we notice absorption in the line core 
and emission in the near line wings (x ~ 4) for ^ = 0° and 10°. This is the characteristic 
self reversal observed in spectral lines formed in extended spherical atmospheres. The self 
reversal decreases gradually as 6 increases, and finally vanishes for large values of 6. Indeed 
for extreme value of ^ = 90°, we observe a pure emission line. 

In Fig. [5]^b) we show line profiles formed in a semi-infinite spherical medium. The 
model parameters are same as in Fig. [5]^a) except for {R, T) = (300, 10^^). The profiles for 
a range of 6' = 0°, 31°, 54°, 72°, 84° are shown. For the core rays {6 = 0°) we see a pure 
absorption line due to thermalization of source function. For other angles, as expected we 
see chromospheric type self-reversed emission lines, formed in the lobe part of the spherical 
medium. 

7. Conclusions 

In this paper we propose a robust method called Preconditioned Bi-Conjugate Gradient 
(Pre-BiCG) method to solve the classical problem of line transfer in spherical media. This 
method belongs to a class of iterative methods based on the projection techniques. We 
briefly present the method, and the computing algorithm. We also present a transpose-free 
variant called the Stabilized Preconditioned Bi-Conjugate Gradient (Pre-BiCG-STAB) 
method which is more advantageous in some of its features. The Pre-BiCG and Pre-BiCG- 
STAB methods are validated in terms of its efficiency and accuracy, by comparing with 
the contemporary iterative methods like Jacobi, GS and SOR. To calculate the benchmark 
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solutions we use spherical shell atmospheres. Few difficult test cases are also presented to 
show that the Pre-BiCG and Pre-BiCG-STAB are efficient numerical methods for spherical 
line transfer. 

L. S. Anusha likes to thank Dr. Han Uitenbroek, Dr. A. Asensio Ramos and Dr. M. 
Sampoorna for useful discussions. 
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Fig. 1. — Geometry of the problem showing the computation of radiation field in a spherically 
symmetric system. The set of core rays and tangent rays are marked. The core is defined 
as a sphere with radius r = 1 in units of the core radius Rcore- The surface is a sphere of 
radius r — R in units of i?core- The rays that intersect the core are called 'core rays', and 
the rest are called 'lobe rays'. No radiation is incident on the outer surface of the sphere 
(outer boundary condition). For all the examples presented in this paper, we use a reflecting 
boundary condition at the z — vertical axis (the mid- line), namely same inner boundary 
conditions are used for both the core and the lobe rays. 
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Fig. 2. — Dependence of the Maximum Relative Change Rc on the iterative progress for 
different methods. Panels (a), (b), and (c) represent models with low, medium and high 
spatial resolution respectively. The model parameters are (n, i?, T, a, e, /3c, By)= (0, 10, 
10^, 10~^, 10~^, 0, 1). The convergence criteria is chosen arbitrarily as (D = 10^^. The 
SOR parameter u =1.5. The figures show clearly that Jacobi method has the smallest 
convergence rate, which progressively increases for GS and SOR methods. Pre-BiCG and 
Pre-BiCG-STAB methods generally have the largest convergence rate compared to the other 
three. 
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Fig. 3. — Behaviour of the true error Tg in a spherically symmetric medium. The model 
parameters are (n, i?, T, e, /3c, B^)— (0, 10, 10^, 10~^, 0, 1). Panel (a) shows the decrement 
of the true error of the Pre-BiCG method for three different spatial grid resolutions. Notice 
the plateau in the true error. Panel (b) shows Te for Pre-BiCG-STAB method. Panel (c) 
shows for different methods, the number of iterations required to reach a constant true error 
2.9 X 10~^. The SOR parameter oj =1.5. The overall behaviour of the curves is same for all 
the methods, although the rates of decrement are different. 
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Fig. 4. — In panel (a) the source function variation with optical depth is shown for a spherical 
media with inverse square opacity variation for two dif f erent values of e. The symbols show 
the benchmark solution read from 



Kunasz fc Hummer 



(jl974l ). which compare well with the 



solution by our method (Pre-BiCG - full lines). The plane parallel solution [R = 1) is shown 
for comparison. Panel (b) shows the effect of power law opacity indices n on the source 
function variation with r. In panel (c) the effects of spherical extension R are shown by 
taking a difficult case of highly scattering, effectively optically thin medium. 



-36- 



o 



(a) n=2 



(b) 7^=10^2; n=2 







-10 



15 



■ /?=lo^ 6=0° ^=1o^ 


9=0° : 


^\ '\ 


6=10° 


^ \ '• \ ' 

\ \ \ 


\ \ \ 7?=1o^ 

\\v 


0=10° 


\\^=1o^ ( 


3 = 90° 




6=96" 





2 4 6 8 10 12 
Frequency x 




2 4 6 8 10 12 
Frequency x 



Fig. 5. — Angular dependence of emergent intensities in highly extended spherical media 
(i? = 10^ and R—lQi^). In panel (a) the profiles for the central ray {9 — 0°), and the lobe 
rays {9 = 10°, 90°) are shown. Panel (b) shows the line profiles formed in a semi-infinite 
spherical atmosphere with R — 300 and T — 10^^. 
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Table 1: The sensitivity of different iterative methods to the convergence criteria uj. Tables 



1(a), 1(b), 1(c) correspond respectively to u;=10 ^, 10 ^, and 10 Number of points per 



decade in the logarithmic r scale is denoted by [A^pts/-D]. The SOR parameter used is 1.5. 

The entries under each method indicate the number of iterations required for convergence. 

(a) w=10~^ 

[A^pts//^] Jacobi GS SOR Pre-BiCG Pre-BiCG-STAB 

5 81 40 24 16 12 

8 136 69 22 19 15 

30 444 230 74 33 23 

(b) (5=10-** 

[A^pts//^] Jacobi GS SOR Pre-BiCG Pre-BiCG-STAB 

5 110 54 30 18 13 

8 186 94 30 22 15 

30 635 325 103 39 30 



(C) £l.= 10 



-10 



[A/^pts//^] Jacobi GS SOR Pre-BiCG Pre-BiCG-STAB 

5 138 68 37 20 14 

8 236 118 40 25 18 

30 827 419 132 45 30 
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A. Construction of A matrix and Preconditioner matrix M 

In Pre-BiCG method, it is essential to compute and store tlie matrix. A brute force 
- fully numerical way of doing this is as follows. Suppose that the dimension of A matrix 
is Ad X Ad, where Ad is the number of depth points. By sending a 5-source function Ad 
times, to a formal solver subroutine. Ad columns of A matrix can be calculated. But this 
takes a large amount of CPU time especially for large values of Ad. 

Instead, there is a semi-analytic way of calculating the A matix. By substituting the 
5-source function in the expression for the intensity on a short-characteristic stencil of 
3-points (MOP in standard notation) we can obtain "recursive relations for intensity matrix 
elements = 1,2,... Ad), which can then be integrated over frequencies and angles to 

get the A matrix. Finally, A = [/ — (1 — e)A]. The diagonal of A is A* and the diagonal of 
A is the preconditioner matrix M. 



Kunasz fc Auerl (119881 ). In radiative 



For plane parallel full-slab problem, this is given in 
transfer problems with spherical symmetry, it is sufficient to compute the solution on a 
quadrant. However this causes a tricky situation, in which we have to define a mid-line 
(see Fig. [T]) on which a non-zero boundary condition = I~ has to be specified. For 
the outgoing rays, the mid-vertical line is the starting grid point for a given ray. Since 

Table 2: Timing efficiency of the iterative methods 

Jacobi GS SOR Pre-BiCG Pre-BiCG-STAB 

CPU time for convergence 7 min 49 sec 4 min 4 sec 1 min 18 sec 27 sec 42 sec 

Overrates in computing 6 sec 6 sec 6 sec 9 sec 6 sec 

Total computing time 7 min 55 sec 4 min 10 sec 1 min 24 sec 36 sec 48 sec 
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the intensity at the starting point is non-zero (J"'" = / ), intensity at any interior point 
depends on the intensity at all the previous points. Recall that 

4(a* > 0) = /k+i(A* > 0) exp(-Ark(/x > 0)) + *k-i5^k-i + ^k^^k + *k+iS'k+i (Al) 

and 

Jk+i(At > 0) = Jk+2(/i > 0) exp(-ATk+i(// > 0)) + ^kS'k + ^k+iS^k+i + *k+2Sfk+2 (A2) 

and so on until we reach the mid-line. It is easy to see from above equations that intensity 
calculation at a short-characteristic stencil MOP is not confined only to the intensity 
on MOP, but also on all previous points, through spatial coupling. This is specific to 
performing radiative transfer on a spherical quadrant. Note that even for the construction 
of a diagonal A, all the elements /y of the intensity matrix has to be computed. We present 
below the recursive relations to compute = 1,2, ... , N^). 

For the incoming rays {jj, < 0) - Reverse sweep 

DO i^l,2,...,Nd 

Consider an arbitrary spatial point i. The delta-source vector is specified as 

S{n)^l, -S(rj) = for i^j. (A3) 
DO ip — 1,2, ... , Np, where A^p is the total number of impact parameters. 

For the inner boundary points, define Ny^ — for the core rays and Ny^ — Nd — {ip—Nc — l) 
for lobe rays. The index Ny represents the total number of points on a given ray of 
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constant impact parameter p. The external boundary condition has to be taken as zero for 



For those rays (with index ip) for which p{ip) < r{i) 

DOj = 2,3,...7V,, 

IF 

(i = i + I) and (j = N^^ or p{ip) = r{j)), which are interior boundary points 
hii^i^^b^^P) = Aj-i('^i'^j-i'^'P)exp(-Arj(/i)) + *u(j,a;,p,/i < 0) + *d(i, < 0) 



This is because at these interior boundary points we assume = Sy, and 8^ = 8^^ = ! 
when j — i + 1. 



constructing integral operators like A. 



Iil{Ti,Ti,X,p) = 0. 



(A4) 



(A5) 



ELSE 



(Non interior boundary points) 



Aj(^i'^j'^'P) = -fij-i(T"i,7j-i,a;,p)exp-(Arj(/i)) + *d(j>,p,/x < 0) 



(A6) 



Elseif j = i 



Ii,i{rurj,x,p) = /ij_i(ri,rj_i,a;,p)exp(-Arj(/i)) + *o(j>,p,/x < 0) 



(AT) 



Elseif j = i - 1 



/i,j(Ti, Tj,x,p) = lij-iin, Tj_i, x,p) exp (-ATj(//)) + x,p,fx< 0) (A8) 
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Else 

if {j<z- 2) 

else 

end if 
End if 
END IF 
END DO 



For those rays for which p{ip) > r{i) 
DO J = 2,3,...7Vvp 



/ij(ri,Tj,a;,p) = (A9) 



Ii,iin,Ti,x,p) = Iij_i(Ti,Tj_i,a;,p)exp(-ATj(/x)) (AlO) 



Ii,j(Ti,Tj,X,p) =0 (All) 



END DO 
END DO 
END DO 

For the outgoing rays (// > 0) - ForwEird sweep 

Let i^Nd,Nd 



S{n) = l, 5(rj)=0 for t^j (A12) 



DO^p= l,2,...,Arp 

For j — Ny^ (Inner boundary point) 
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Ii,iin,rj,x,p) (/X > 0) = Iij(ri,rj,a;,p) (/x < 0) 



(A13) 



For non boundary points 



For those rays (with index ip) for which p{ip) < r{i) 



B0j = N,^-l,N,^-2,...l 

Ii,j(ri, rj,x,p) = /i,j+i(Ti, Tj+i, x,p) exp (-ATj(Ai)) + *d(j, x,p,ii> 0) (A14) 

Elseif j = i 

= -fij+i('^i'Ti+i'^'P)exp(-ATj(/x)) + *o(j>,P,/^ > 0) (A15) 

Elseif j = i — 1 

Ii,i(Ti,Tj,x,p) = /ij+i(Ti,rj+i,a;,p)exp(-ATj(/x)) + *u(j,a;,p,/x > 0) (A16) 



Else 



Ti,x,p) = Ii,j+i(ri, Tj+i, exp (-Arj(//)) 



(A17) 



End if 



END DO 




Iij(Ti,Tj,X,p) = 



(A18) 



-43- 



END DO 
END DO 
END DO 

The algorithm given above saves a great deal of computing time by cutting down the 
number of calls to the formal solver -2 instead of A^d -the first call to store the * and At 
at all depth points, and the second call to compute lij 



